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Abstract. We consider elongated condensates that cross twice the speed of sound. In 
^fi ' the absence of periodic boundary conditions, the phonon spectrum possesses a discrete 

and finite set of complex frequency modes that induce a laser effect. This effect consti- 
£2 \ tutes a dynamical instability and is due to the fact that the supersonic region acts as a 

tJjT) 1 resonant cavity. We numerically compute the complex frequencies and density-density 

correlation function. We obtain patterns with very specific signatures. In terms of 
the gravitational analogy, the flows we consider correspond to a pair of black hole and 
r-; ■ white hole horizons, and the laser effect can be conceived as a self-amplified Hawking 

^^| radiation. This is verified by comparing the outgoing flux at early time with the stan- 

dard black hole radiation. 
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1. Introduction 

In 1981, Unruh suggested [1] that the analogue of black hole radiation could be observed 
in a quantum fluid that crosses the speed of sound, thereby forming a sonic horizon. 
Since then, several types of fluids and several types of flows were considered [2]. Very 
recently, a near-stationary supersonic flow engendering a pair of horizons (a black hole 
one followed by a white hole one) was realized in a Bose-Einstein condensate [3]. In 
such a background, because of the anomalous dispersion of Bogoliubov phonons, one 
expects to get a kind of laser effect [4]-[6] due to a self-amplification of the Hawking 
radiation. 

The first aim of this paper is to provide the theoretical basis of this effect starting 
from the Bogoliubov-de Gennes equation, i.e. without making use of the gravitational 
analogy. In this we apply to flows containing two horizons the treatment of [7] that 
was only applied to a single (black hole or white hole) sonic horizon. Following [6] we 
then establish that the laser effect is governed by a discrete set of complex frequency 
modes that correspond to the resonant modes of the cavity formed by the region 
bordered by the two horizons. In a classical description, this discrete set governs 
the dynamical instability of the flow. The second aim is to compute the spectrum 



of the complex frequency modes by theoretical and numerical methods. The excellent 
agreement of the results validates both the concepts and the semi-classical methods used 
in the theoretical approach. Finally we consider two observables: the mean number 
of emitted phonons, and the two-point function of the phonon field that governs the 
density-density correlation pattern [8, 9]. At sufficiently late time, both observables 
are governed by a single mode, the most unstable one. In this late time regime, their 
behaviours are identical to those one would obtain using a classical description of the 
density perturbations. At early time instead, when starting from vacuum configurations, 
correspondence with the quantum phonon flux emitted by a black hole horizon [7] is 
established. 

The propagation of phonons in flows containing two horizons has already been 
considered. However, toroidal configurations with periodic boundary conditions were 
generally used [10, 11]. In that case, the spectrum is very complicated because it results 
from a combination of two effects: the discreteness of the wave vectors defined on the 
torus interferes with that associated with modes that are trapped in the supersonic 
region. As a result, not only the analysis is difficult, but the relationships with the 
Hawking effect and the black hole laser effect [4]-[6] are hard to draw. On the contrary, 
when dealing with continuous wave vectors, the analysis of the complex frequency mode 
is simpler, and the relationship with Hawking radiation is easily made. 

We have organized this paper as follows. In section 2 we present the Bogoliubov-de 
Gennes equation in a way that is suitable for our aims. In section 3 we determine the 
spectrum of asymptotically bound modes, and explain why complex frequency modes 
appear when the speed of sound is twice crossed. We then compute the eigen-frequencies 
using semi-classical methods. In section 4 we numerically solve the Bogoliubov-de 
Gennes equation and compute the spectrum. We then evaluate the mean occupation 
number and compare it with the spectrum one would obtain if only the black hole 
horizon were present. We finally evaluate the correlation pattern of the density-density 
two-point function. 

2. Framework 

2.1. Black-hole-white-hole geometries 

Using the analogy [1, 2] between the wave equation governing sound propagation in 
a moving fluid and that governing light propagation in a curved space-time, one can 
associate an acoustic geometry to the fluid flow. In 1+1 dimensions, these geometries 
are of the form [4] 

ds 2 = -c 2 dt 2 + (dx - vdtf, (1) 

where v is the flow velocity and c the speed of sound. In this language, a black-hole- 
white-hole geometry is obtained when v crosses twice c. Indeed, a sonic horizon is 
present whenever v crosses c. Assuming that the fluid flows from right to left (v < 0), 
the location x H of a horizon is given by c(xn) + u(#h) — 0. In our case, we set the 



location of the white hole horizon and that of the black hole respectively to iw = — L 
and xb = L. These divide the x-axis in three regions, which we call I on the left of the 
white horizon (x < —L), II between the two horizons (— L < x < L) and III on the 
right of the black horizon (x > L). The flow is supersonic in the internal region II and 
subsonic otherwise. The opposite case where the velocity is subsonic inside does not 
give rise to a laser effect, and will be considered elsewhere. 

In the present work, our analysis is restricted to the following stationary profiles: 

Ky?\x + L\ 



c(x) + v (x) = ch-D sign(x 2 — L 2 ) tanh ' n 



tanh 1/n 



Kb \x 



cnD 



,(2) 



cnD 

which generalize to two horizons those used in [7]. ch is the sound speed at the horizons, 
2L is the distance between them, D determines the size of the near-horizon regions where 
the metric is not flat, n controls the sharpness of the transition to the flat regions, and 
/t\v and «b control the surface gravities. Indeed, the surface gravity is defined by [2] 

1 d(c 2 - v 2 ) 



il< 



dx 



(3) 

X=Xn 

which yields -ch«w for the white hole and ch^b for the black hole. The metric (1) is 
then completely fixed by introducing an extra parameter q: 

c(x) = ca + (l-q)[c(x)+v(x)], 
v(x) = -c H + q[c(x) + v(x)], 
which specifies how c + v is shared between c and v. When restricting to the cases where 
c{x) > 0, v(x) < 0, D > 0, (5) 

the range of D and q is limited. The allowed values are graphically illustrated by the 
shaded area in figure 1, left panel. Three velocity profiles corresponding to different 
couples (q, D) are plotted in figure 1, right panel. 

2.2. Density perturbations in Bose Einstein condensates 

We present the main steps leading to the equations for linear density perturbations in 
Bose-Einstein condensates [12], having in mind cases where the condensate crosses the 
speed of sound. We follow [7], where more details can be found. 

At low temperature, the quantum properties of a gas of weakly interacting atoms 
are efficiently described by a second quantized field satisfying the commutation relation 

[%x),*t( i , x ')] =5 3 (x _ x / )) (6) 

and by its Hamiltonian 



H= [d 3 x\ — V x ^ t V x ^ + y^ t ^ + -^ t ^ t ^^l 
J [2m 2 J 



(7) 



where m is the atom mass, V the external potential and g the effective coupling. The 
last two quantities can depend on both t and x. When a significant fraction of the 
atoms condense, it is meaningful to expand \1/ in a c- number function ty , describing the 
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Figure 1. Left panel: allowed values for q and D giving rise to a black- hole-white-hole 
geometry. Right panel: Velocity profiles c(x)/ch (black lines) and |w(x)|/ch (red lines) 
for three different values of (q,D): q = 0.8, D = 1 (solid lines), q = — 1, D = 0.3 
(dashed lines) and q = 2, D = 0.1 (dotted lines). These values are also reported in the 
left plot. All profiles with % = kb = iz. 



condensed part, plus a field operator <f>, describing (relative) density perturbations over 
the condensate 

# = *o(l + 0)- (8) 

In what follows, we consider elongated condensates, which means that the 
transverse excitations have sufficiently high energies that they are not excited. In this 
case, \l/o an d 4> are effectively one- dimensional fields. For simplicity, we also assume 
that the condensate is infinitely long, in order to avoid discussing the discreteness of the 
longitudinal wave number k. The discrete case is indeed more complicated, and treated 
in details in [11]. Finally, we assume that the condensate is stationary, which means 
that \l/o is °f the form 

qr (t, x) = e-'^ t/h x v/p (x)e ie ° (x) , (9) 

where \x is the chemical potential, po(x) the mean density of condensed atoms and 

(10) 



h 

v(x) = —d x 6 (x) 
m 



is their mean velocity In this case, the Gross-Pitaevskii equation 

h 2 



ma* 



t^o 



■£fi + V + g* 



\Pr 



reduces to 



/' 



-mv 



h 2 d 2 xy fp~ 



V + gpo, d x (vp ) = 0, 



2 2m po 

where the second equation is the continuity equation for a stationary flow. 



;n) 



(12) 
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At the linear level, obeys the Bogoliubov-de Gennes equation. In the present 
case, using (8) and (12), this equation reads 

ihd t (j) = [T p - ivhd x + mc 2 ] + mc 2 ^, (13) 

where we introduced the x-dependent speed of sound 

c 2 (x) = 9{X)P0{X \ (14) 

m 

and where 

T >=-L^ %pod * = -L vd 'l d ' (15) 

reduces to the usual kinetic operator when the background condensate is homogeneous. 
The second expression is only valid in stationary flows. Equation (13) tells us that the 
perturbation is coupled to the condensate only through the functions v and c. The 
disappearance, from (13), of the potential V, the quantum potential — ^- x v po , and the 
coupling g constitutes an essential step towards the notion of the acoustic metric of (1); 
see [7] for more details. 

Using the canonical commutation relation 

[4>(t,x),ft(t,x')} = J—Six-x 1 ), (16) 

Po{x) 

which follows from equations (6) and (8), (13) is the Heisenberg equation 

m t 4> = [4>,H cS \, (17) 

engendered by the Hermitian Hamiltonian 

H cS = fdx ^ j> [T p - ivhd x ] + [{T p + ivhd x )ft]$ + g Po (<j> ]2 + ft + 20^) } , (18) 

which is obtained by expanding (7) to second order in 0, 0t; see appendix A. 

3. Theoretical analysis 

3.1. Spectrum of bound modes and quantization 

In stationary condensates, the spectrum of (18) can be characterized using very general 
properties. These are first, the fact that the eigenmodes are asymptotically bound 
(to guarantee that they have a well-defined norm since the spatial domain of is 
infinite), second, the Hermiticity of (18), i.e. its self-adjointness with respect to the 
scalar product defining the eigen-modes norm, and third, that this scalar product is not 
positive definite, as it is the case here, see (B.15), and for bosonic fields, see e.g. [13]. 

When these conditions are met, the spectrum of (18) contains a continuous set of 
real frequency modes labelled by to and a discrete index a, plus, possibly, a discrete and 
finite set of complex frequency modes that appear in pairs with complex conjugated 
frequencies. (If the scalar product were positive definite, as is the case for fermions, the 
spectrum of (18) would be purely real.) Throughout this paper, complex frequencies 



shall be written as X a = u a + iV a , where a is a positive integer which labels the discrete 
set of pairs, and where u a and T a are both real and positive. The other cases are reached 
by complex conjugation and/or multiplication by —1. The discrete index a describes 
the subset of modes with the same real frequency. As explained below, in one dimension, 
it contains two or four modes depending on whether the flow is subsonic or supersonic. 
As a result, when propagating in a stationary condensate, 4>, obeying (16) and (17), 
can always be expanded (see appendix B) as 



4>(t,x) 



du 



Z^ 



K iwi €(*K 



+ e 



+iuit 



(vSC*))^] 



Yl [e~' lXat Ux)ba + e- iA °Va(x)c fl + e +ix ^(r] a (x)rbl + e +iA »*(Ca(x))*ct] , (19) 



where the modes satisfy the normalization relations 



dx Po [(rMniix) - Mx)ytffi(x)\ = 5 aa ,5(u - u/) 
dxp [(£ a (x))*ilJa>(x) - (r) a (x))*( a > (x)] = i5 Xa \,. 



(20) 
(21) 



Correspondingly, the operators a, b, c satisfy the commutation relations 

= 5 aa '8(u 

iS 



1 ui a u}' \ 



UJ 



(22) 
Kc[,} = i6 XaKl . (23) 

All the other scalar products, and all other commutators, vanish. In particular, one 
has [b a , b a ] = [c a , c* a ,] = 0. Hence b a and c a are not destruction operators. In fact, 
the operators b a and c a form pairs [6] that characterize complex, i.e. with two degrees 
of freedom, harmonic oscillators that are unstable, and for which therefore there is 
no notion of quanta. (Concomitantly, one verifies that the norm Jdx p [|£a| 2 — |?7a| 2 ] 
vanishes.) It should be noticed that the eigenfrequencies associated with b a and c a are 
complex conjugated from each other (respectively A a and A*), something guaranteed by 
the Hermiticity of (18). 

In terms of the above eigenmodes and operators, (18) reads 



H, 



cff 






a a ?a a , 



blc a 



dx Po (x)\vZ(x)\' 



dx po(x)( a (x)(rj a (x)y 



h.c. 



(24) 



The first line contains the usual sum over harmonic oscillators of real frequency and 
the c-number term accounting for the depletion, while the second line is due to the 
complex frequency modes. It should be noticed that the unusual form of the second 
c-number term follows from the unusual scalar product of (21). This Hamiltonian is not 
bounded from below due to the b\c a terms. Hence the vacuum can no longer be defined 
as the ground state of H, as one might have expected since one is dealing with unstable 
oscillators. This gives rise to some ambiguity when choosing the initial 'vacuum' state 
(see section 3.4). 
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Eqs. (19) and (24) can be rewritten in a more familiar form by decomposing each 
couple (b a , c a ) into two couples of destruction/creation operators d a+ , d a+ and d a _, d„_, 
as shown in appendix B. The field becomes 



fax) = [dcoJ2 [e^rMdZ + e^i^x))*^] 

J a 

+ y^ yj> a+ {t,x)d a+ + 4> a _(t,x)d a - + (cp a+ (t,x))*di + + (tp a -(t,x))*dl_ 

a 

and the Hamiltonian (24) 

H cS = JduJ2^ (V«2 - J dx Po (x)\^(xy 2 



(25) 



a 






d 1 d 



dxp (x) (\cp a+ (t, x)\ 2 - \cp a _(t,x)\ 2 ) 



dl + rfl_ + / dxp (x) ((f) a+ (t,x)<p a -(t,x)) - h.c. 



(26) 



We conclude with two remarks. Firstly, if one can describe the discrete set using 
the d a operators, there is a price to pay as the associated modes 4> a+ , <p a +, 4>a-, Pa- are 
not frequency eigenmodes. Hence their time dependence cannot be factorized. Secondly, 
in usual circumstances, the discrete set of complex frequency modes is empty, as is the 
case in homogeneous and in near-homogeneous condensates. Nevertheless, when the 
condensate crosses twice the speed of sound, the discrete set is (generally) nonempty. 
To explain why, we first need to study the spectrum in homogeneous subsonic and 
supersonic flows, and then understand how to paste the corresponding modes when the 
flow crosses the speed of sound. When solving this, the dispersive properties of become 
essential, as we now recall, and as first understood in [4]. 

3.2. Mode analysis 

Inserting the field expansion (19) in (13), one obtains a system of two c-number equations 

[h(X + ivd x ) -T p - mc 2 ] A = mc 2 ip x , 

[— h(\ + ivd x ) —T p — mc 2 ~\ tp\ = mc 2 (f)x- 

The couple (<f>\,(p\) can be formed by either the real frequency modes (^2(x), ip"(x)), 
or the complex frequency modes (£ a (x) , r] a (x)) or (i(; a (x), Ca{%)) an d, accordingly, A can 
be real or complex. Eliminating ip\ from the above system, one obtains [7] 



[h(X + ivd x ) + Tp] - [-h(\ + ivd x ) + Tp 



h 2 vd x —d a 

v 



0. 



(28) 



When the background quantities are independent of x, Fourier modes 4>\ oc exp(ik\x) 
are solutions of (28), provided k\ is a (possibly complex) root of the dispersion relation 

C H 4 fc 4 



(A - vk) 2 = c 2 k 2 + 



A 2 



n 2 (k), 



(29) 





cak/n ^it<^ cnk/K, 



Figure 2. Graphical solution of the dispersion relation (29) for subsonic flow (left 
panel) and supersonic flow (right panel). Solid lines: ±Q(k)/K. Left panel, dashed 
line: (uj — vk)/n. Right panel, dashed, dotted, dotdashed lines: w — vk respectively for 
to < CJ m ax! w = w max , uj > uj max . For uj < uj max , one clearly sees that, in supersonic 
flows, two extra (real) roots exist in the left lower quadrant. The most negative one is 
called, in the text, ku and the other one kb , so that fci, ' < ku < 0. 



where Q is the frequency in the comoving frame, and where ch is a typical value of the 
sound speed, we take to be that at the horizons, see section 2.1. We also introduced 

A = — f~, (30) 

which gives the characteristic dispersive scale. When sending A — > oo, (29) becomes 
the relativistic equation Q 2 = c 2 k 2 since the quartic term drops out. Similarly, (28) 
becomes the Euler equation governing sound waves. Instead, when keeping this term, 
the dispersion relation (29) possesses four roots, some of which can be complex. As we 
shall progressively see, the two extra roots are at the origin of the laser effect. 

3.2.1. Modes in homogeneous condensates. In subsonic flows, \v\ < c, for real X = u, 
two roots are real, and describe the right and left moving solutions. The other two are 
complex, conjugated to each other, and correspond to modes that are asymptotically 
growing or decaying, say to left. Only the first two should be used in (19), as the last 
two are not asymptotically bound. 

In supersonic flows, the situation is quite different. For real X = u>, there exists a 
critical frequency w max such that the four roots of (29) are [7] 

• uj < w max : all real; 

• uj > w max : two real and two complex ones, as in subsonic flows. 

See figure 2 for a graphical solution of the dispersion relation (29). 

When u is real, the normalization of the (bound) modes can be easily worked out, 
as it is, up to a trivial factor, independent of the constant velocity v (because of Galilean 



10 
invariance). Let us write 

€(x)e~^ = uZ, <{x)e~^ = t£, (31) 

V27TPO V 2 7rpo 

where a spans over four values if \v\ > c and u < u ma , x and over two values elsewhere. 

Using the mode normalization (20) 

(l<| 2 -K| 2 )^ = l, (32) 

from the mode equation (27) and the dispersion relation, one has 

Dks< = C (33) 

where 

1 



D, 



mc 2 



h k h k 2 

h\ I c 2 k 2 + - — t mc 

Am z 2m 



(34) 



is independent of v. Equations (32) and (33) fix u^ and w" except for a common phase. 
Taking both w" and t>" real, one obtains 



(35) 



In supersonic flows, when < u < w max , the two modes associated with the extra real 
roots of (29) have negative norm. However, for — CJ max < uj < 0, the corresponding 
modes have a positive norm. So, the field is expanded as$ 





/0fc° 1 e ifc 2 x 


/<%£ ,D fc2 e ifc S- 


V <9w A _ D 2 a V2VTP0 



+e +iwi | [M*^ + M*a* f ] + %- - w) X) 0-L«-L | J , (37) 

where 6(z) is the Heaviside function. 

In subsonic flows one obtains the same expression, but without the last term in 
each bracket since the extra real roots are absent. 

3.2.2. Modes in black-hole-white-hole geometries. To obtain the spectrum in 
condensates that cross twice the speed of sound, i.e. in flows defining the geometry 
of (2), we need to take into account the nontrivial propagation in the non- homogeneous 
regions localized near the two horizons, and the fact that each subsonic region is now 

| In terms of the doublets W denned in appendix B one has 



oo 



W = I dco [W^al + WZal + W£3tf + W%a v J] + I dw V \W®a®, + W^a^ . (36) 



o Jo 



»=1,2 



11 

defined on the half line and no longer on the full real axis. Then, we need to identify 
the number of globally defined independent modes. 

It should be noted that this set does not depend on the particular form of the 
mode equation, but mainly resides on the behaviour in the internal region II and in 
the asymptotic flat external regions I and III described in section 2.1. Therefore, the 
result is the same when studying a phonon field in a BEC, as in the present paper, or a 
generic scalar field with superluminal dispersion relation as in [6]. As explained in this 
reference, in geometries described by (2), there is a continuous double set (right-going 
and left-going waves) of real frequency modes, and the discrete set of complex frequency 
modes is generally not empty. Here we give a more mathematical argument yielding the 
same result. 

In flows (2), two asymptotically flat and subsonic regions play the most important 
role in defining the set of modes, for both A real and A complex. For A = u real, there 
are three bounded modes: the two usual propagating modes associated with the real 
roots of (29) plus the mode that asymptotically decays. This mode should now be taken 
into account as it is bound in the subsonic region where it is defined. We will call 
the two sets of three modes 0^ ! and 0^ m , respectively for the left (I) and right (III) 
asymptotic region, where the superscript % takes three values: u,v,CJ to characterize 
the asymptotic right-going w-mode, the left-going f-mode and the decaying CJ-mode, 
named after the paper of Corley and Jacobson [14]. In the internal region II, we have 
four modes that we call <^ n with j — 1,2,3, 4. Since this region is bound, the four roots 
of (29) should be taken into account. 

Therefore, any globally defined solution of frequency u can be expanded in two 
forms by referring to its left, or right, asymptotic behaviour: 



8,1 

U) Vol ) 



i 



m 



Moreover, each partial wave 0^ T or 0^ m can be written in term of the 0i n : 

' (39) 

3 

Putting together the above equations, one obtains 

This system of four equations in six unknowns (the coefficients R^ and L % J) has a two- 
dimensional set of solutions, corresponding to two linearly independent modes. For 
instance, one can choose as independent solutions, 0^ ,m and <p v u ] in , the two right-going 
and left-going in-modes, defined as the solution of (40) with respectively L u = 1, R" = 
and L u — 0,R V — 1. Similarly, one can chose the out-modes (^ out and (j)^j out defined 
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respectively by Li" = 0, R u = 1 and L v = 1, R u = 0. Using the scalar product of (20) to 
normalize these modes, the two sets are related by [6] 



rp ruout I D J.V,OUt 
v,out i d iu,out 



(41) 



Unitarity imposes |TJ 2 + |/L| 2 = 1 = (TLl 2 + l^] 2 , and i^T* + T^* = 0. 

Let us now move to the complex frequency case. Because all expressions are analytic 
in A, the above analysis applies as such when replacing u by A = u + iT. The novel 
aspects only come through the condition that the modes be asymptotically bound. 
Indeed, T > implies that fc", the asymptotic wave number of the u real root of (29) 
(more precisely of its analytical continuation in T), acquires a positive imaginary part. 
Thus the mode diverges at x — > — oo, unless one puts L u = 0. Similarly, on the right 
side, the v mode diverges at x — > +oo, and imposing that it is bound requires R v = 0. 
However, these two conditions imply that the system (40) has only the trivial solution 
R l = L % = 0, except when its determinant vanishes. This condition defines an equation 
for the complex frequency A that has a finite set of solutions: {A a , a = 1,2,..., iV}. 

3.3. The semi-classical treatment of [6] 

Having established the condition that defines the complex frequencies A a , we now turn 
to the question of calculating them. Since the above algebraic method does not depend 
on the specific form of the mode equation, the semi-classical treatment of [6] applies 
to the present case. Therefore, we just summarize the results without going into the 
details. One should nevertheless pay attention to the fact that one is dealing with modes 
doublets (<f), ip) and to the normalization of these modes. 

When the background is not homogeneous but v varies slowly with respect to the 
wavelength of the perturbation, the exact solutions are well approximated by their WKB 
approximation, which can be directly inferred from (35) 



<%{*) 



</£(*) 



dk£(x) 




1 


exp 


[i f x dx'k%(x 


')] 


duj 


D 


- n 2 

U kZ(.x 
k%(x) 


) 
exp 


^2irp (x) 
[i J x dx'k"(ji 




'dk-(x) 


0] 



(42) 



9U a/ 1 " D kl*) ^^ {X) 



The ^-dependent wave numbers k"(x) are real roots of the dispersion relation (29) in a 
stationary inhomogeneous flow characterized by v(x) and c(x): 

(oo-v(x)k) 2 = c(x) 2 k 2 + ^\. (43) 

4m 2 

Ignoring the quartic term, one obtains the dispersion relation of a massless field 
propagating in the acoustic metric (1). Equations are valid in the three regions I, 
II, III but not close to the horizons where the gradients of v and c cannot be neglected. 
The theoretical analysis can then be greatly simplified when taking into account 
the fact that the mixing between u and v modes is usually negligible [7, 15]. (However 
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the numerical analysis of the next section does not rely on this simplifying assumption.) 
In this approximation, the laser effect is completely due to the u-modes (for flows to 
the left v < 0). Considering the u in-mode in the internal region II, it can be written 
as a superposition of three WKB waves: 

C in = A,C + CVS* + CVS*, (44) 

where Lp_ w are the two negative frequency y9 u -modes, corresponding to the two extra 
roots hi , , fci, being the most negative one, see figure 2. 

In [16] it was shown that the propagation from one horizon to the other is efficiently 
described by a unitary scattering of a two-component vector, whose first component 
represents the w-wave 0^ and the other represents the trapped wave, described either 
by (if_l)* or (ip_l)*. The 5" matrix is decomposed in four elementary matrices 

S=U 4 U 3 U 2 U U (45) 

representing (from 1 to 4) the scattering at the white hole horizon, the propagation 
between the white and the black hole horizon of 0^ and (cp_i)*, the scattering at the 
black hole horizon and the propagation back to the white hole horizon of (ip_^)*. These 
matrices are 



(46) 



where 

pL pFluj 

1,2, (47) 

are the Hamilton-Jacobi actions governing the phase of the WKB modes, and L u and 
Ru are the two turning points. Moreover, by unitarity, the parameters in U\ and U 3 
must satisfy |aj 2 = |<SJ 2 , I7J 2 = \%\ 2 , |« w | 2 (l - \z w \ 2 ) = |7^| 2 (1 - \w LU \ 2 ). 

For real frequency modes, the condition that the trapped mode be single valued 
imposes [6] 
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where U is real by unitarity. From this matricial equation, one obtains 

S21 e i^ = _|lil^||2 ; (49) 



1 — O22 >->22 1 — ^22 

and therefore the coefficients of (44) are 

A u = a u (l + z u b u ), BW = a(i£ + b u ), B® = b u . (50) 

For complex frequency with T > 0, imposing that the incoming u branch vanishes, 
one obtains 
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which implies 

& = S 12l S 22 = 1. (52) 

Equations (49), (50) and (52) are the central equations we shall use to analyse 
the results of the numerical integration. Even though these equations have been 
obtained through a semi-classical reasoning, their validity goes beyond that of the 
WKB approximation, as shall be established by the remarkable agreement between 
the relations one can derive from them and the numerical results. 

The main prediction concerns the relation between A a , solutions of (52), and the 
behaviour of the coefficients of (50). Indeed, A a solution of S 22 = 1, is a pole of 
Buj — bu = £21/(1 — S22). Therefore, when T a is small enough (this is always true 
in our numerical situations), the pole is close to the real axis and the behaviour of Bh 
for real frequencies u is dominated by the contribution of the pole: 

»(2) _ $21 lF ° «(2) (z.o\ 

u ~l~S 22 ~ u -( Ua + T a ) B »" {b6) 



from which 



?(2)|2 



&?» ,„ _^, + P lggl'- («) 



Thus, \Bw I should be well described by a sum of Lorentzians characterized by the 
complex frequencies A . In figure 3 we plot \Bu \ 2 obtained from numerical simulations 
and the corresponding fitted sum of Lorentzians. The agreement is excellent. 

We conclude this analysis by showing that a semi-classical treatment furnishes 
approximate expressions for the complex frequencies A . Assuming |£ w | 2 , \w 2 \ and 
l^w^wl are niuch smaller than 1, one can expand the condition S 22 = 1 in powers of 
these quantities. To zeroth order, one has A a = u a real, with u a solution of the Bohr- 
Sommerfeld condition [6] 

A„ 
S { -i a ~ s -l a ~ arg(aa, a 7^J = / dx [-k£}(x) + k {2) a (x)] - arg(a w 7 w ) = 2vm B s, (55) 

where u-qq, = a — 1,2, . . . ,N. A semi-classical treatment gives arg(a^7 w ) = — n. A more 
refined estimate [17] gives 

arg(a^7j « -7T - / I — ) - / I — ), 

\«b/ \k-wJ (56) 

f{x) = arg T E (ix) — x log(x) + x + — , 

where Te is Euler's Gamma function. 

To first order in |^| 2 , \w 2 |, |2 w u> w | and 5\ a = 5u a + iT a , S 22 = 1 gives 

2T*Y a = (\z u f + K| + 2|^ <l c os^ a ) = \S 12 M\ 2 , (57) 

T u a SuJ " = -\zw a w* Ua \ sin^ a , (58) 

# a = S: a + sU a + M g Z "« W }° a "\ (59) 



Ua duj L 



c( 2 ) c(!) 1 ~ ~ 



(60) 
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Figure 3. We represent |f?£, | 2 as a function of w/cj max for two values of L, 
Lk/ch = 6 (upper panel), Lk/ch = 25 (lower panel), and for % = kb = k, q = 0.5, 
D = 0.33, n = 1, A/k = 2, u max /K ~ 0.195. Green points: numerical simulation; 
red lines: fitted series of Lorentzians. Online movie eigenfrequencies.gif (available 
from stacks. iop.org/NJP/12/095015/mmedia): Evolution of \Bu | 2 from Lk/ch = 15- 
25. 



where T^ is the time for the trapped mode to make a full bounce. The phase d a 
modulates the frequency imaginary part T a and the first order correction to its real part 
5u a , as shall be seen in section 4.2.2. 



3.4- Observables: Phonon fluxes and density-density correlation patterns 

Having identified the set of eigenmodes, we now consider observable quantities. We 
shall study both the phonon flux emitted by the black-hole-white-hole system, and the 
non-local density-density correlation pattern, as they illustrate very different aspects of 
the laser effect. 

In what follows, we work in quantum settings because we shall assume that the 
state at the formation of the supersonic region is vacuum. In classical settings, the 
dynamical instability (the laser effect) would be present only if the initial density- 
perturbation possesses a non-vanishing overlap with some complex frequency mode, 
see section V.A. in [6]. For a perfectly adiabatic formation of the supersonic flow, the 
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amplitudes associated with these modes would be zero. However, any small deviation 
from adiabaticity will trigger the instability, thereby engendering a behaviour very 
similar to that derived in quantum settings. In fact, as we shall see, the main difference 
between the quantum and the classical treatment is at early times because in the 
quantum vacuum all complex frequency modes contribute to the observables through 
the spontaneous excitation. 

To compute observables, we first need to choose the quantum state. Because of the 
instability, the energy (24) is unbounded from below, and there is no clear definition 
of the vacuum. Nevertheless, if the formation of the supersonic region at time to = 
is adiabatic, and if the temperature condensate is low enough, i.e. much smaller than 
k, the Heisenberg state of the system can be approximated by the state annihilated, at 
that initial time, by the destruction operators a", d a+ , d a - of (25). § Since the operators 
d a +, d a - are not stationary, the expectation values will not be stationary either, and 
will instead depend on the lapse of time since the formation at to = 0. In what follows 
all expectation values will be computed in this state. 

The two observables we shall study are related to the density fluctuation p\ = p—po- 
Given the definition of in (8), fi\ is given by 

h = Po(0 + f ) = poX, (61) 

where x — + $ is Hermitian. Its expansion in terms of operators a, d is 

X{t,x) = J<kuJ2[ e ~ i " t xZ( x K + h - c -]+ I] [xa+{t,x)d a+ + X a-{t, x)4_ + h.c.] , (62) 

a a 

where 

x2(aOMS(aO + ¥£(*), 

Xa,± {t,x) = (p a ,± (t,x) + <fa,± (t, x) . (63) 

In the vacuum defined at to = 0, the two-point function is 
(0\ X (t,x)x(t',x')\0) = [dtu^-'^-^xZWix^x'))* 

J a 

+ J2(xa + (t,x)(Xa+(t',x')y + Xa-(t,x)( X a-(t',x')y). (64) 

a 

Because of the complex frequency modes, it is not a function oit — t' only. This differs 
from the cases of a single black hole (or white hole) horizon, where the two-point function 
is stationary in vacuum [7]. 

To extract more physical information, it is convenient to return to the frequency 
eigenmodes appearing in (19). Defining 

&a(x) = £a(x) + Va(x), U a (x) = lp a (x) + (a(x), (65) 

§ We applied the standard reasoning to both the a" and the d operators because r a , the imaginary 
part of the complex frequency A a , is smaller than a tenth of u a = RcA a . 
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the two-point function becomes 



(0\x(t,x)x(t',x')\0) = fdooJ^e-^-^xZWixZW))* 

J a 

} r(«+*')e- iw -(*-*') - (x)(o- (x'))* + e- r(t+t ' ) e- ,w -( t - t ' ) i/ B (x)(i/ a (x')) 

- iA ^V a (x)K(x')r - <r iX W-^v a (x){a a {x')y] . (66) 



+E Re 



+i y^ Re U 



The first line gives the (vacuum) contribution of the real frequency modes. Since these 
are only elastically scattered, see (41), they stay in their ground state and contribute 
neither to the emitted fluxes nor to the pattern of density-density correlations. The third 
line contains a mixed contribution of the growing a a (x) and decaying modes v a {x)- It 
does not contribute to the fluxes at any time because cr a (x) vanishes in region I, and 
v a {x) does it in III, see the paragraph after (41). Moreover, since it only depends on t—t', 
it does not significantly contribute to the correlation pattern at late time. Therefore, 
both fluxes and correlation patterns are governed by the growing mode contribution, 
the first term of the second line in a a (x)(a a (x'))*. 

It is clear that at late times, i.e. t — to > l/max(r a ), all observables are dominated 
by the mode with the largest T a . For instance, the equal time density-density correlation 
function is asymptotically given by 

(OlpifoaOpifoaOlO) ~ Po(x)po(x') x e 2r "'Re [a a (x)(a a (x'))*} . (67) 

The real part of the frequency u a drops out and, as a result, the locus of the maxima 
does not propagate with time. This function is studied for various modes in section 4.4. 
In case the initial state is not vacuum but a thermal state, the above two-point function 
would be multiplied by 2n a + l, where n a is the mean occupation number, which depends 
in a nontrivial way both on the temperature and on the blue shift effect due to the 
horizons (see equation (46) in [7]). Similar considerations apply to the quantity that 
is considered below. Had we worked in classical settings, we would have obtained a 
coherent wave whose behaviour in x and t is the same as that of (0\pi(t,x)pi(t',x')\0) 
at fixed x',t', see equation (54) in [6]. The main difference arises from the fact that the 
classical wave amplitude is fixed by initial conditions. 

It is also interesting to study the onset of the instability for earlier time, i.e. 
t — to < l/max(r a ). To this end, we consider the following quantity: 

P bWh (u;,T)= f dU I dte- 1 ^-^ (<d\ X {t,x) X {t\x)\Q), 
Jo Jo 

v r Q T, / N | 2 4|sin[(cj-q; a -zr a )T/2]| 2 
~ E " 6 MX)1 (.-^ + (T a )i • (68) 

It governs the probability that a density fluctuation of frequency u be observed during 
the interval [0,T] at some fixed location x in the subsonic region III, see [6]. 

When r o — > 0, when the discrete set labeled a becomes continuous, and when T 
is sufficiently large, Pbh-wh(^,T) behaves as in the Golden Rule: Pbh-wh(^,T) oc Tn w , 
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i.e. the lapse of time T times n^, the mean occupation number of frequency u. There 
can be pre-factors that depend on the strength of the coupling, the amplitude of the 
mode 0^ at x, if one deals with a derivative coupling. 

It is more interesting to study Pbh-wh(w, T) as a function of T for a given black- 
hole-white-hole geometry, and to compare it with Pbh(w, T), the corresponding quantity 
computed in the case where only the black hole horizon would be present. When 
T x max(r a ) <C 1, the two observables behave in a very similar manner, even when 
only a few (say 10) complex frequency modes exist in the black- hole- white-hole case. 
This shall be seen by a numerical analysis in section 4.3, but can be also understood 
from the properties of the growing modes cr a (x), see the central paragraph of section 
V.B.5 in [6]. 

For classical settings, the laser effect is present only when the initial density- 
perturbation profile is characterized by a non-vanishing amplitude of some complex 
frequency mode. Since for adiabatic formation of the horizons all the amplitudes 
associated with those modes are zero, no laser effect appears in this case. However, 
even a small deviation from perfect adiabaticity causes the presence of instabilities. 
The main difference in the quantum framework is therefore the possibility to make the 
system unstable through the spontaneous excitation of complex frequency modes. 

4. Numerical results 

4-1. The method 

In a few words, we describe how we proceeded to solve numerically (27) in the flows 
described by (2) and (4). As in the case of a single black hole or white hole horizon, 
when working with a fixed frequency, the main task is to avoid the growing mode 
contaminating the oscillatory modes. The way to get rid of this difficulty is by 
integrating from the subsonic region towards the horizon into the supersonic region 
where there is no growing mode [14]. 

Basically we solved separately the mode equation from region I to region II, and 
from region III to region II, using a code adapted from [7]. In each case, as in 
that reference, we integrated the equations for initial conditions describing the three 
acceptable modes in regions I and III: the right-moving w-mode, the left-moving -u-mode 
and the decaying mode. Then for each of them we computed the amplitudes of the four 
oscillatory modes in the supersonic region II. The two globally defined modes are built 
using the procedure described in section 3.2.2. To obtain the scattering from I to III, 
we eliminated the amplitudes in region II. 

To obtain the correlation patterns of figures 11-13 we solved the mode equation 
with the corresponding complex frequency X a that we had formerly computed. 



19 

4-2. The discrete set of complex frequencies 

Two strategies can be used to determine the complex frequencies A a = u a +iT a , solutions 
of (52). They can be obtained either by solving directly the linear system (40), or by 
determining the center and the width of the Lorentzians in \B U | 2 of (54). We use the 
latter to localize them and the former to refine the results and study how they depend 
on the various parameters of the system. 

To give a first idea, in figure 3 we represent \Bl | 2 as a function of u/u maiX (green 
points) for two different values of L, the distance between the horizons. A sum of 
Lorentzians functions (red line) is fitted to the data obtained from the numerical analysis 
(see (54)). The quality of the fit confirms the correctness of the theoretical analysis 
of section 3.3 and of [6]. In the following sections, we study the dependence of A a 
on both the parameters describing the geometry (see the velocity profile (2)) and the 
dispersive scale A of (29). The dependence on n is not reported here, since it induces 
no significant change. || Before proceeding we make some comment about the properties 
and the validity of the Bohr-Sommerfeld condition (55). 

4-2.1. The validity of the semi- classical approximation. The action appearing in (55) 
contains the difference of two wave vectors: kl — kl . These have the same sign 
(negative for positive w), as they belong to the same w-branch of the dispersion relation 
(see figure 2), but they have opposite group velocity: kl describes a right-going mode 

(2) 

with respect to the lab, whereas k u a left-moving one. This peculiarity give rise to 
an unusual phenomenon. In the usual case, eigenmodes with small (large) frequency 
correspond to small (large) Bohr-Sommerfeld numbers u-qs, i-e to modes with few 
(many) nodes. In the present case instead, a large nes corresponds to low frequency 
modes and vice versa. This can be understood from (55). Because the wave vectors 
appeared subtracted, a small nes implies ku close to k u and this happens for ui close 
to cj max (see figure 2). On the contrary, a large 7Zbs implies a large difference between 
ku and k w , that is, u close to 0. 

This has an unusual consequence. The Bohr-Sommerfeld condition is more reliable 
when the action is large, that is for ubs 3> 1- In the present case this happens for small 
u/u max . On the other hand the WKB approximation is expected to fail for small u, as 
can be verified by the fact that / of (56) cannot be neglected when u/k < 1. Both these 
expectations are confirmed in figure 4 where we compare the numerical results with the 
predictions obtained with the standard WKB approximation (green lines) and with the 
improved method (red lines), that is when using (56). Firstly, the agreement between 
the Bohr-Sommerfeld condition and the numerical results is worse for to close to cj max 
and gets better when nes increases. Secondly, at low frequency, while the quality of the 

The role of n of (2) is to govern the smoothness of the transition from the near-horizon region to the 
flat asymptotic ones. A smaller n corresponds to smoother transition, while a larger n corresponds to 
steeper transition and, consequently, to a larger almost-flat region between the horizons. The usefulness 
of n is essentially technical: it allows one to control the numerical analysis when L, the distance between 
the horizons, is comparable to D, the width of the transition between regions I— III of section 2.1. 
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Figure 4. \Bu | 2 as a function of w/w max for % = Kb = K, q = 0.5, D = 0.5, n = 2, 
Lk/cu = 10, A/k = 8, w max /« ra 1.408. Points: numerical simulation; green dashed 
spikes: standard WKB approximation; red solid spikes: improved treatment of (56). 



standard WKB prediction (green lines) becomes worse, the improved method continues 
to work very well, thereby establishing its validity. Further studies about the agreement 
of the improved method and numerical results are in preparation [17]. 



4-2.2. The birth of modes as a function of the distance between the horizons. In figure 5, 
left panel, we plot uj a = Re(A a ) for the entire spectrum of complex frequency modes as a 
function of the half-distance between the horizons L, all other quantities being fixed. On 
the right panel, T a = Im(A a ) is plotted as a function of L, for ?t,bs=4 and for the same 
values of the other parameters. As conjectured in [6], when L grows, all u a (L) values 
increase, and when there is enough 'room', a new eigenmode appears with u a ~ 0. As 
a result, the eigenfrequencies become denser close to w ma x, because they neither cross it 
nor disappear. 

It is interesting to consider more closely the birth of new modes near u = 0. A 
rather difficult question to answer is the following: when extra mode appear, does 
the imaginary part of the frequency T a vanish, as found in figure 1 of [18], or not? 
This question is physically relevant in that it governs the continuous character of the 
stability of the system: if T a vanishes, it implies that the birth of new modes leads to a 
continuous behaviour in L (as one can expect on the general ground that the properties 
of the eigenmodes, solutions of (13), continuously depend on the parameters appearing 
in that equation). From figure 5, which is based on the roots of the determinant in (40), 
no definite conclusion can be drawn since there is a loss of numerical precision for 
u) — > 0. A more reliable method consists in studying the behaviour of \Bw | 2 on 
the real frequency axis. We created a gif-animation (eigenfrequencies.gif (available 
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Figure 5. Left panel: The set of aj a /w max as a function of L, from nBs=l (upper 
curve) to 71bs = 12 (lower curve) for the same fixed parameters as in figure 3. At 
kL/ch = 25, the values of w a /oj max correspond to the peaks of the lower panel of 
that figure from right (tt-bs = 1) to left (tt-bs = 12). The broad peak of figure 3 near 
u/w mal[ = 0.05 cannot be reproduced with this method because of numerical errors. 
For similar reasons, each curve w a (L) ends for low u>. Right panel: T a /n for ubs = 4 
as a function of L for the same fixed parameters. The occurrence of zeros is due to the 
fact that kb — kw, as discussed in the text. 



from stacks. iop.org/NJP/12/095015/mmedia)) of \Bu | 2 as a function of a; for increasing 
values of L to illustrate the appearance of the new eigenfrequencies. From this it is quite 
clear that when a new eigenmode appears, both the real part and the imaginary part 
of A a vanish, as was also found in the appendix of [19]. 

Another important feature confirming the theoretical analysis of section 3.3 is the 
presence of superposed oscillations on the overall trend of u a . In fact, the latter is given 
by the solution of the Bohr-Sommerfeld equation (55), whereas the oscillations are due 
to the deviations governed by the sine in 5u a of (58). For the imaginary part instead, 
there is no zeroth-order contribution, and the dominant contribution to the right plot 
of figure 5 is given in (57). 



4-2.3. Asymmetric velocity profiles. Using asymmetric velocity profiles, namely 
different surface gravities n w and Kb, does not lead to substantial modifications of 
the spectrum. The only observable change concerns T a . When kw = ^b, T a vanishes 
for some particular combination of parameters, as can be seen in figure 5, right panel. 
Equation (57) can indeed be rewritten as 



1 



2T„ 



(koj - KJ) 2 + 4 l^a<J 



COS 



•&« 



(69) 



which shows that when 



z. 



\w„ 



T n vanishes when the cosine does, and this is 



"U> a I I "'Wa I 1 

because the scatterings at the black and the white horizons destructively interfere with 
each other. When the two surface gravities are different, (\zu a \ — I^J) is non-zero 
and T a can no longer vanish (see figure 6). The non-zero offset is significant only for 
sufficiently large values of cu/k, which requires that u m3iX /n be large enough, since all 
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Figure 6. Left panel: T a /n-B as a function of Lkb/ch for Kw = k b (solid line), 
kw/ k b = 0.5 (dotted line) and kw/«b = 2 (dashed line), for w max /«B = 0.974, 
L> = 0.33, n=l,q = 0.5, A/k b = 10. 

Using the standard expressions z u = e~ n0J / Kw and w u = e~ 7rw / K s ) which 
furnish reliable estimates [7], the ratio between the offset and the amplitude of the 
oscillation is 

*2 

= smh 2 (iru a (K B — k w )/2). (70) 



w u 



J) S 



4 l^a<J 



Even though this function increases in u when K\v 7^ «b, the consequences of this 
are never important because when the ratio is large, the corresponding mode will not 
significantly contribute to the instability since r a <C 1, as can be seen from (69). 



4-2.4- The u^v mixing and the parameter q. The mixing between the w-modes, which 
are right-going in the lab in subsonic flows, and the f-modes is controlled by the 
transmission coefficient T u of (41). In figure 7, \T U \ 2 is plotted as a function of q 
for a fixed frequency oj/k = 0.100. The transmission coefficient is almost 1 between 
q = 0.25 and q = 0.75. Thus, for values of q in this range the u-v mixing is negligible, 
as was noticed in [7] for a single black (or white) hole. In this regime the approximation 
discussed in section 3.3 is valid. There is, however, an important difference with respect 
to the single black hole case. In black-hole-white-hole geometries, the transmission 
coefficient deviates from zero when approaching a resonance. When working with a 
fixed u, we found two sharp spikes figure 7 for q = 0.1248 and q = 0.56023 which 
exactly correspond to two unstable modes. Moreover, the width of the spikes is almost 
equal to the corresponding T a . We also looked for complex eigenfrequency associated 
with the other two local and broad minima of |T| 2 , at q — 0.3363 and q = 0.6285, and 
we found two eigenmodes with respectively u} a /n = 0.090 and uj a /n = 0.103, i.e. in 
the neighbourhood of the chosen frequency u/n = 0.100. These considerations show 
that close to resonances one cannot completely neglect the u-v mixing. Nevertheless, as 
shown in figure 4, the improved treatment of the Bohr-Sommerfeld equation produces 
very good estimates for the values of the real part u a of the complex frequencies. 
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Figure 7. Transmission coefficient |T U | 2 as a function of q, at constant ui/k = 0.100, 
and for the other parameters as in the bottom plot of figure 3. The sharp minima 
at q = 0.1248 and q = 0.56023 correspond to complex frequencies at oj a j 'k = 0.100, 
whereas the broad ones at q = 0.3363 and q = 0.6285 correspond to frequencies 
respectively at uj a / k = 0.090 and u) a / k = 0.103. 



4-2.5. The role of the maximal frequency w max . In [7, 15], it was shown that the 
deviations, due to dispersion and w.r.t. the standard Planckian distribution, of the 
spectrum emitted by a single BH, or WH, are mainly governed by w max . However, the 
latter depends both on the UV scale A and the velocity profile (see section 3.2.2). The 
relationship takes the form 

u max = Af(D,q). (71) 

Hence the same value of w max can be reached from very different cases, and yet it was 
found that the fluxes are hardly sensitive to this. In the present case however, this 
insensitivity is lost because the number of resonances directly depends on A. This can 
be understood by studying (55), and is manifest in figure 8, where \B^ | 2 is plotted as 



a function of oj/k for two different values of A. 



4-3. Growth of the asymptotic phonon fluxes 

In figure 9, we represent the relevant quantities that govern the properties of the 
'probability' function introduced in (68). On the left upper plot, we present the growth 
rates T a and the corresponding times T^ of (60) for the 13 complex frequency modes that 
exist in the condensate flow of figure 3, lower plot. On the right upper plot, we show the 
Golden Rule transition rate associated with (68) in that black-hole-white-hole geometry 
(dots), and the rate one would obtain if the white hole were not present. In that case, 
the quantity plotted is proportional to w x n u , where n^ is the mean occupation number 



in the black hole geometry as computed in [7]. It is given by n u 
see Us of (46), where {w^l 2 is well approximated by 



w u 



V(i 



w u 



w„ 



-^/"(l-uVuw) 1 / 2 . 



(72) 
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Figure 8. \Bu | 2 as a function of w/w max for different values of D and A giving rise 
to the same w max /ft = 0.19: black solid line for D = 0.33, A/k = 2, and red dashed 
line for D = 0.7, A/k = 0.695. The other parameters are those of the bottom plot of 
figure 3. 
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Figure 9. Upper left panel: bouncing time Tj of (60) as a function of u>/u> max in 
unit of k" 1 (points) and corresponding values of T a (boxes), when Lk/cjj = 25. The 
condensate parameters are as in figure 3, lower plot. Upper right panel, solid line: 
transition rate associated with (68) due to the radiation emitted by the sole black hole 
as a function of w/w max (solid line); corresponding quantity due to the discrete set of 
complex frequency modes in the black-hole-white-hole geometry. The oscillations are 
due to the cosine in (57). Lower panels: P{uj,T) of (68) as a function of w/w max for 
an isolated black hole (solid line) and for a black- hole-white-hole pair (dashed line), 
respectively at time T = 30k -1 (left plot) and T = 200k -1 (right plot). 
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The pre- factor of u is due to the normalization of the x modes of (63). In a black-hole- 
white-hole, the quantity that corresponds to u x n^ is u a x 2T a T^ . In both cases, these 
quantities vanish for uj > u> max . 

In the lower plots, we show (68), at two different times, and for both the black- 
hole-white-hole, and the isolated black hole flows. As expected, the exponential growth 
of the laser effect does not show up for times T smaller than the inverse of the maximal 
T a (here T ~ 110/k). Moreover, the discreteness of the spectrum is not visible either 
at early times. To be resolved, it requires times larger than 2nAu a , where Au a is the 
frequency gap between neighboring Bohr-Sommerfeld frequencies u a . It is here of the 
order of Au a ~ cj max /10 ~ 0.02k. 

On the contrary, for times of the order of l/max(r a ) or greater, both the exponential 
growth of the laser effect and the discreteness of the spectrum show up. These allow 
one to distinguish the phonon flux emitted by the black-hole-white-hole pair from that 
emitted by the sole black hole that grows linearly in T for all values of u. This linear 
growth can be observed by comparing the continuous lines of the right upper and lower 
plots, and constitutes a numerical validation of the Golden Rule! 

4-4- Spatial properties of complex frequency modes and correlation patterns 

The density-density correlation function can be calculated following the procedure 
described in section 3.4. As outlined there, only the largest-r a mode significantly 
contributes to the pattern at late time. Nevertheless, to appreciate the variety of cases, 
it is worth studying the modes and the corresponding correlation patterns also for lower 

r. 

As a typical example of a mode with a large T, we consider the fourth BS mode in 
a configuration with fivecomplex eigenfrequencies. In figure 10 (upper panel) the real 
part of this mode is represented, while the supplementary movie laser_l.gif (available 
from stacks. iop.org/NJP/12/095015/mmedia) gives its evolution in time. As a second 
example, we choose a very different situation with q = 0.7 and a dispersion relation 
with higher A. In this case, the spectrum of complex frequencies is large. In figure 10 
two modes are plotted: one with tt-bs = 14 and a moderate value of T (central panel), 
and one with ubs — 2 and a very small value of T (lower panel). When comparing the 
modes and their correlation patterns, some features remain the same, whereas others 
significantly differ. 

In the first example, since T is high, the laser effect is strong as can be seen from the 
movie: in the right asymptotic region there is an exponentially growing right-going flux. 
This growth in time can also be inferred from the spatial decrease of the mode amplitude 
on the right of the black hole, because these are linearly related, see equation (54) of 
[6]. Considering the correlation pattern of this mode in figure 11, one sees that the 
dominant signals come from the inside region (the central square of size ~ 22) where 
the amplitude of the mode is highest, and from the horizontal and vertical bands which 
describe the correlations between the escaping flux and the trapped mode. 
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Figure 10. Real part of the growing mode £, a (x) of (19) for respectively high, moderate 
and low values of I\ Red lines: white and black horizons. Upper panel: fourth 
cigenmodc with lu^/k = 0.098, T^/k = 0.0056, with Lk/ch = 10.95 and the other 
parameters as those of figure 3. Central panel: eigenmode with ubs = 14, wu/k = 0.46 
and Th/k = 0.004. w max /ft = 2.53, q = 0.7, D = 0.7, n = 2, % = kb = «, 
A/k = 8, Lk/ch = 5. Lower panel: eigenmode in the same background with tibs = 2, 
UJ2/K = 2.48, T2/ k = 2. x 10~ 6 . The importance of T can be estimated from the relative 
mode amplitude in the right outside region versus that inside in the trapped region. 
Online movies showing the evolution of ^ a (^) as a function of t, respectively laser_l.gif 
(from t = -400/k to 400/k), laser_2.gif (from t = -200/k to 200/k), laser_3.gif (from 
t = -400/k to 400/k) (available from stacks.iop.org/NJP/12/095015/mmedia). 
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Figure 11. Density-density correlation function (67) divided by e 2r<1 * of the mode 
with a high T represented in the upper plot of figure 10 with x in units of Cn/n. We have 
divided the expression by e 2Tat in order to suppress the dependence on time. Since the 
correlations bewtcen different regions have very different amplitudes, we were obliged 
to use different scales to reveal them. The central square from — 1 1 till 11 is the trapped 
region. It contains the strongest correlations since the mode amplitude is highest in 
it. The two (symmetrical) bands of width ~ 22 in the right upper plot describe the 
correlations between the trapped region and the external region on the right of the 
black hole, while the other two bands describe those with the left external region. The 
correlations in the lower left panel are those on the left of the white hole. They are 
much weaker and of greater wave length, in agreement with the mode properties on 
the left side of the upper plot of figure 10. 
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Figure 12. Density- density correlation function (67) divided by e 2Tat for the mode 
with a moderate T represented in the central plot of figure 10. Similarities with the 
former figure are manifest, and concern both the spatial properties of the pattern and 
the amplitudes of the correlations. The main difference concerns the high-frequency 
modulations of the inside-inside and inside-outside correlations which are due to the 
fast oscillations of the mode amplitude in the trapped region, which arc visible in the 
central plot of figure 10, and which are due to fact that nes is high: ubs = 14. 



From the mode itself and from its correlation pattern in figure 11, one can see 
that the u-v mixing is very low, as expected from the choice of q = 0.5. The left- 
going mode coming out of the white hole is about 2.5 orders of magnitudes smaller 
than the right-going one that comes out of the black hole, as can be seen from the 
ratio of amplitudes between the top left part and top right part in figure 11. To 
understand how the modes propagate, it is appropriate to consider the top left part. 
This region represents the correlations between u (right-going) and v (left-going) modes. 
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Figure 13. Density- density correlation function (67) divided by e 2Tat for the mode 
with a low T represented in the lowest panel of figure 10. The bottom legend 
corresponds to the central region of the plot. One clearly sees the two nodes in the 
trapped region. One also finds them in the inside-outside correlations, both on the left 
and on the right of the trapped region. Finally, one notices that in the present case 
the u— u correlations on the right upper square are weaker than the v — v correlations 
on the lower left panel. This follows from the high u-v coupling between the trapped 
mode, which is a u-mode, and the outside w-modc. 



The slope of the highest /lowest correlation lines gives the ratio between the velocities 
of the two modes. When the modes form a continuous set, the two velocities are the 
group velocities. However, in the present case, when considering a single mode, the 
slope of these lines corresponds instead to the ratio between the phase velocities of u- 
and -y-modes. Similarly, the pattern in the top right part are simply due to the nodes 
of the right moving mode evaluated at a given time, see equation (61) in [6]. 

We now consider the correlation pattern of figure 12, corresponding to a mode 
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with a moderate value of T and ?t,bs very large. The overall properties of the pattern 
are basically the same as those of figure 11. The main differences come from the high 
wave- vector content of the trapped mode, which produces short distance features in the 
central square and in the bands. Outside these regions the patterns are very similar. 

In the last pattern of figure 13, corresponding to a mode with a very small value 
of T and 77-bs = 2, we get a different picture. The pattern is basically concentrated in 
the central square since the amplitude of the mode outside is related to its amplitude 
inside by the square root of T, as can be seen from (57). In the present case, the u-v 
mixing is so high that the pattern on the left of the white hole is higher (by a factor of 
about 5) than on the right of the black hole. One also notices that in the central panel 
(ubs — 2) there is only one node. Moreover, the beat-like shape shows that this wave is 
the result of the superposition of two waves with very similar and relatively high wave 
vectors. 

4-5. The Technion experiment 

Recently a black-hole-white- hole flow has been experimentally realized [3]. Our 
numerical code can describe this configuration, up to some limitation. In fact, the 
velocity profile (figure 3(b) of [3]) is rather irregular while our code works accurately 
only for velocity profiles as in figure 1, right panel, with an almost-flat internal region. 
Nevertheless, we shall proceed in order to obtain at least an estimation on the number of 
the unstable modes and their growing rate. As far as we know, this is the first estimate 
of this quantity, which is very relevant from an experimental standpoint. 
The main parameters of the experimental setting are 

c H = 7.3 x 10" 4 m/s, 

Kb = Kw/2 = k = 213.7 s~ , 

v (x = 0) = -3.6 x 1(T 3 m/s = -4.9 c H , 

c( x = 0) = 3.5 x 10" 4 m/s = 0.48 c H , (73) 

c(x ->• oo) = 8 x 10^ 4 m/s =1.1 c H , (74) 

2L = 10" 5 m = 1.5c h /k, 

A/k = 6. 

To reproduce the above setting we perform two simulations with fixed D = 1. 
In figure 14, \Bl | 2 is plotted as a function of w/a; max for q = 0.9 (left panel) and 
q = 0.48 (right panel). These two values are obtained by using, respectively, c(x = 0) 
and c{x — > oo) equal to their experimental values in (73) and (74). ^f The small distance 
between the two horizon implies that few trapped modes are present in this situation. 
The values of T a and u a of the complex eigenfrequencies are reported in both cases 
in table 1. The dominant contribution comes from the eigenfrequency with the largest 

% We used these two values because, in our analysis, it is not possible to fix independently c(x = 0) 
and c(x — > oo). We expect that the actual value of the complex frequencies would lie in between the 
two sets we obtain. 
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Figure 14. |f?£, | 2 as a function of w/w max for q = 0.9, u> max / k ps 3.4 (left panel) and 
g = 0.48, w max /«; ~ 2.7 (right panel). The values of the other parameters are given 
in the text. Green points: numerical simulation; red lines: fitted series of Lorentzians. 
oj a and r o from the fit are reported in table 1. 



Table 1. lu ci and T a for the complex eigenfrequencies in the systems described in the 
caption of figure 14. 



q = 0.9 



q = 0.48 



w /w n 



uj a /n 



T a /n cj a /w r] 



LO a /K 



r a /« 



2 x 10~ 9 7 x 10~ 9 0.076 2 x 10~ 8 5 x 10~ 8 0.23 

0.12 0.41 0.012 0.32 0.86 0.022 

0.47 1.58 0.005 0.68 1.80 0.004 

0.75 2.54 0.0007 0.92 2.43 0.0005 

0.93 3.17 0.0004 



r a . In the two cases this corresponds to an instability time scale, respectively r « 0.06 s 
and r rs 0.02 s. Even if the two simulations are generated with very different values 
of q, the two results agree in order of magnitude. Note that the estimated instability 
time scale is more than twice the time scale for which the horizons are maintained in 
the experimental configuration (« 0.008 s). In order to see the laser effect, it would be 
needed to maintain the systems for a longer time, or to increase k by a factor of 10, 
without changing the adimensional parameters of the system such as the ratio Lk/cu- 

5. Conclusions 



In this paper, we analysed the spectrum of phonons in a BEC when the stationary 
flow of the condensate crosses twice the speed of sound. Our analysis is based on 
the Bogoliubov-de Gennes equation (13). Hence, even though the analogy with light 
propagation in a pair of black and white horizons is manifest, none of our results rely 
on this gravitational analogy. 

In the limit where the condensate can be considered as infinite, i.e. when no periodic 
boundary condition is introduced, the spectrum of bound modes contains real frequency 
modes that are only elastically scattered, see (41), plus a discrete and finite set of pairs 
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of complex frequency modes (19). These modes can be seen as the resonances of the 
cavity bordered by the two sonic horizons. They lead to dynamical instabilities because 
the scattering through each sonic horizon is anomalous, in that it mixes positive and 
negative norm modes, see (44). Then, using semi-classical techniques, we compute the 
real and the imaginary part of the complex frequencies. The real part u a obeys a 
Bohr-Sommerfeld condition that introduces the discreteness of the spectrum, whereas 
the imaginary part T a is related to the norm of the scattering coefficients across the 
horizons, see (57). In the case of toroidal geometry, the analysis would instead be 
complicated by the discreteness of the wave vectors. The instability appears only when 
one bj a approximately matches one of the frequencies associated with that discrete set 
of wave vectors. This explains the presence of the instability bands found in [10, 11]. 

In section 4, we numerically solved the Bogoliubov-de Gennes equation for real 
and complex frequencies. By comparing the numerical properties with those derived 
using the Bohr-Sommerfeld and (44), we validated the use of semi-classical methods, 
see figures 3-5. In section 4.3 we numerically studied the growth of number of phonons 
emitted by the black-hole-white-hole system. In spite of the discrete character of the 
trapped modes, at early time, this quantity behaves very similarly to the flux that 
the black hole horizon would emit in the absence of the white hole horizon. Instead, 
for larger time both the instability and the discreteness show up. Finally, we studied 
the equal time correlation pattern of density-density fluctuations associated with three 
different complex frequency modes, respectively with high, moderate and low value of the 
imaginary part T a . In all cases, the instability shows up most clearly in the supersonic 
region where the amplitude of the trapped mode is the highest one. The correlations 
between the trapped mode and the emitted modes display very specific properties, see 
figure 11 and then next two ones. 

Finally, we studied (within some limitations) the experimental situation realized in 
June 2009 in the Technion [3]. The results are summarized in figure 14 and show that 
few unstable modes are found, and that the instability time scale is about ten times 
larger than the life time of the condensate. This means that an increase by a factor of 
ten of the surface gravity (without changing the adimensional parameters of the system 
such as the ratio Lk/ch) would lead to comparable time scales. The laser effect could 
then be observable. 
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Appendix A. The Hamiltonian of density perturbations 

We expand H of (7) in powers of <fi defined in (8) to the second order: 

h 2 



H = jdx% 



;m ^ + fp * r 



H 1 = dx %ft 
H 2 = [dxpolft 



~dl + v + m 



^o + h.c. 



ft 2 (9 2 *n 

T - ivhd x - 2-± + V + 2gp 

2m W 



9 



Pofvp 



where T p is defined in (15). Expanding the Heisenberg equation of motion for ^/, 



iW t ^(t,x) = [if(t, x), H], 



one obtains 



m*o = * [4>,H 1 ] 



ihdf 



Ho 



ihxj) 






The former one gives the Gross-Pit aevskii equation (11). Defining 

h 2 d 2 j> Q 



5H 



dxpo 



2m \1> 
equation (A. 6) can be rewritten as ihd t < 

Hctt = H 2 -5H = / dx p { (f) j [T p - ivkd x + gp 



V + gpo 4> j 4>, 
\(/>,H cS ], where 
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P0(X\2 



(A.l) 
(A.2) 
(A3) 

(A.4) 

(A.5) 
(A.6) 

(A.7) 
(A.8) 



Remarkably, it is independent of external potential V and of c^^o/^o- When 
symmetrizing it with respect to <p and <jy , one obtains (18) which is manifestly Hermitian. 

Appendix B. The quantization procedure 

To understand how to expand the complex field operator in modes and 
creation/destruction operators, it is convenient to define a two-component field [16] 

(B.9) 

(B.10) 
(B.H) 

(B.12) 
(B.13) 



igpo<72, 



w = 

Then (13) becomes 

ihd t W = BW, 

B = (T p + gp )<7 3 -ivd x 

where Oi are the Pauli matrices 

1\ /0 

i oj' a2= {i 

Since the field W is invariant under the conjugation operation defined by 
S = a x S\ 



<7\ 



&3 



1 
-1 
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the structure of W must be 

W = J2( w nan + W n at l ), (B.14) 

n 

where W n are doublets of C-functions, and J2 n denotes the summation over a (possibly 
continuous) complete sets of modes. 

Using <9j = —d x , Tj = T p and the properties of Pauli matrices, one verifies that 
the scalar product 



(W 1 \W 2 ) = J dx p (x) W*(t, x)a 3 W 2 (t, x) (B.15) 

is conserved under time evolution when Wi are solution of (B.10), since 

B* T a 3 = a 3 B. (B.16) 

For later convenience, we state here some properties of the scalar product that can 
be verified using the anticommutation relation of Pauli matrices and (16): 

(Wt\W 2 ) = (W 2 |Wi)*, (B.17) 

{W 1 \W 2 ) = -(W 2 \W 1 ) = -(WilWb)*, (B.18) 

WW T - {WW T ) T = —^t5{x - x')ia 2 , (B.19) 

Po(x) 

[(WtlW), (W 2 \W)] = -{W X \W 2 ) = {W 2 \W X ). (B.20) 

Assuming that the W n have a non-zero norm, we define an orthonormal basis: 

(W n \W m ) = - (W n \W m ) = 5 nm , (B.21) 

(W n \W m ) = 0, (B.22) 

where the Kronecker 5 is replaced by a ^-distribution in the case of a continuous set of 
modes. Using (B.14) and the above identities, one obtains 

[a n ,al] = [{W n \W),-{W m \W)} = (W n \W m ) = 5 nm , (B.23) 

which shows that a n and a^ are in fact destruction and creation operators. 

When the condensate is stationary, one can work with eigenmodes W® of frequency 
A, where the index a describes the set of modes with the same frequency. By definition, 
one has 

BW% = H\W?. (B.24) 

The conservation of the scalar product gives 

= d t (WZ\W« f ') = -i(A* - X)(WZ\WZ'). (B.25) 

From this it is clear that only real frequency modes can be normalized as in (B.21). 
However, in the presence of dynamical instabilities, complex frequency eigenmodes are 
present. We shall assume that these modes form a discrete and finite set of pairs of 
modes with conjugated frequencies that we call {A a = u a + fT a , a = 1, 2, ...iV}. In place 
of (B.21), we choose the following pseudo-normalization for each pair, 

(W Xa \W x *,)=i5 aa ,, (B.26) 
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and the other scalar product must vanish because of (B.25). In fact, the Hermiticity of 
H implies that if A a is an eigenfrequency, then A* is an eigenfrequency too [16]. We call 
V\ a and Z\ a the doublets corresponding to A a and A*, respectively. 
Summarizing, the field W can be expanded as 

W = JdcoJ2 [WSK + W^atf] +J2 [Vak + Z a c a + W a + Z a c\] , (B.27) 

a a 

where all A a have a positive imaginary part T a . The modes satisfy the following 
normalization rules, 

(W£|W#) = - (W5\W3) = 8{u - u/)<W, 
(V a \Z a ,) = (V a \Z a ,) = -(Z a ,\V a ) = -(Z a ,\V a ) = i<W, 
and all the other scalar product vanish. Using 

b a = i(Z a \W), 

h{ = i(z a \w), 

C a = -i{V a \W), 

c\= -i{V a \W) 
and (B.20), one obtains 

[S a ,4] = [i(Z a \W),-i(V a ,\W)} = ~{Z a \V al ) = i5 aa ,, 
which is the commutation relation of complex unstable oscillators [6]. Decomposing the 
doublets as 



(B.28) 
(B.29) 

(B.30) 
(B.31) 
(B.32) 
(B.33) 

(B.34) 



W, 



, -kjt l<l%(x)\ y = e _ iAat / Ca(x) A z = e _ iA * t ( i) a (x) 



(B.35) 



the field (f> can be expanded as (19), and equations (B.28) and (B.29) give (20) and (21). 
Using (19), H c f[ is a sum of three terms, containing respectively only real frequency 
modes, only complex frequency modes and mixing real and complex frequencies: 

H e s = H T + H c + H mix . (B.36) 

Using the commutation relations and the scalar products of section 3.1, we obtain 

h 



H v = — dudu' y^ / dx 



, Po{x) 



y + ^e+^C^re-^^R^ 



+ {oo — uj )e 



>\-A-iOJt(j,ai 



(x))*e +i - \<fo{p)Ya*a a J + (u/ - u^^^x)^ VM*) W 



po(x)e 



-(a;' + a;)e"-VS(^)e +iwi (^(^))*C«2' 
5 fdoodoo' J> + cS)(WZ\W2)&$&i 

aa' 

— — / dudu' 2_,( u + u')8(ui — u')5 aa / / dx 

aa' 

+ \ Jdcodu' £>(w«| wtfja^ + \ Jdcodoj' £ 

aa' aa' 

/ do; 2_] hu 



-i(u)'-u)tf( a 



((^(x)T^(x)) 



cj(WS\W$)1Cc& 






dx Po {x)\^{x)f 



(B.37) 
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Using the same techniques, it is easy to show that H m { x = and 



H c = ihJ2{K blca+l 



fe+ / dxp (x)( a (x)(r) a (x))* 



h.c. 



(B.38) 



Putting everything together, one obtains (24). 

Note that 6's and c's are not destruction operators [see (23)]. Nevertheless, they 
can be combined to give couples of creation/destruction operators 

k + ict 



3 - C w ja + Ka 



da- = 



W- u a ~ 10 a 



V2 



(B.39) 
(B.40) 



which in fact satisfy the commutation relations 

[4+,4+] = A-, 4'-] = S aa', 

and all the other commutators vanish. Equation (B.27) becomes 

W = fdco J2 [W<X + W°a a J] +]T [\V a J a+ + W a .d a . + W a+ d\ + + W a .d\_] , (B.41) 

a a 

where 

V a - iZ a 



W a+ = e"^ 



Wa _ = e -»- Va - 



(B.42) 



V2 ' a ~~ V2 ' 

whose normalization is 

(W a+ \W a , + )=(W a _\W a ,_) =-{W a+ \W a , + ) =-(W a _\W a ,_) =5 aa ,, (B.43) 

and the other scalar products vanish. Notice that W a ± are no longer frequency 
eigenmodes. Decomposing W a+ and W a - as 

\V 2 a+(i,a;)/ \(p a -{t,X)J 

the field expansion (19) becomes (25), and the Hamiltonian (24) becomes (26) when the 
(arbitrary) phases 9± are chosen to be 6± = 0. 
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